% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Setting
% BV_Target = {'LB'};                                                        
BV_Target = {'UB'};                                                        
% BV_Target = {'LBUB'};

%% Load Data
FileName = ['AllTradingDate19962021_Monthly.mat'];
load(FileName, ...
     'Target_AllDate');             
clear FileName

Date_Monthly = Target_AllDate;                                
clear Target_AllDate 

%% Optimal Risk Preferences under Rank-Dependent Expected Utility
AllCR_Target = 1;
OPT_TAB = nan(length(AllCR_Target), 6, size(Date_Monthly, 1));             
OPT_OP_LB = cell(length(AllCR_Target), size(Date_Monthly, 1));             
OPT_OP_UB = cell(length(AllCR_Target), size(Date_Monthly, 1));             
EU_OP_LB = cell(1, size(Date_Monthly, 1));                                 
EU_OP_UB = cell(1, size(Date_Monthly, 1));                                 
for d = 1:size(Date_Monthly, 1)
    Target_Date = Date_Monthly(d, 1);
    Target_TTM = Date_Monthly(d, 3);
    
    % Load Data
    FileName = ['Summary_OP_Grid_' num2str(Target_Date)];
    load(FileName, ...
         'Record', 'Record_OP_LB', 'Record_OP_UB', ...
         'OP_Bid', 'OP_Ask', 'IV_Bid', 'IV_Ask', ...
         'K', 'KS', 'CPFlag', 'S0', 'S0_ADJ', 'TTM', 'RF', 'DY');
    clear FileName

    Type = cell(size(K));                                                  
    Type(CPFlag==1) = {'Call'};   
    Type(CPFlag==2) = {'Put'};    

    % Expected Utility
    Index = find((Record(:, 1)==1) & ...
                 (Record(:, 2)==1));
    OP_LB = Record_OP_LB(Index, :)';   
    OP_LB = max(OP_LB, 0);                                                 
    EU_OP_LB(1, d) = cellstr(char(num2str(OP_LB')));                                     
    OP_UB = Record_OP_UB(Index, :)'; 
    OP_UB = max(OP_UB, 0);                                                
    EU_OP_UB(1, d) = cellstr(char(num2str(OP_UB')));        
    clear Index OP_LB OP_UB  
    
    % Rank-Dependent Expected Utility
    for r = 1:length(AllCR_Target)
        CR_Target = AllCR_Target(r);
        NUM = round(CR_Target * length(K));                           
        CR_Real = NUM / length(K);                                    
        clear NUM
        
        % Based on Lower Bound
        if strcmp(BV_Target, 'LB')==1
            Index_Real = find(Record(:, end - 2)==CR_Real);                
            if length(Index_Real) > 0
                Index_EU = find((Record(Index_Real, 1)==1) & ...
                                (Record(Index_Real, 2)==1));     
                if length(Index_EU) > 0
                    Index_Real = Index_Real(Index_EU);                           
                    OPT_TAB(r, :, d) = [CR_Target Record(Index_Real, 1:2) CR_Real nan nan];
                    OPT_OP_LB(r, d) = EU_OP_LB(1, d);
                    OPT_OP_UB(r, d) =  EU_OP_UB(1, d);
                else
                    Record_CR_Real = nan(size(Index_Real));                
                    Record_DEV_Real = nan(size(Index_Real));               
                    for c = 1:length(Index_Real)            
                        OP_LB = Record_OP_LB(Index_Real(c), :)';   
                        OP_LB = max(OP_LB, 0);                             

                        Index_Ask_W_L = logical(OP_LB <= OP_Ask);
                        Index_Ask_W_H = logical(OP_LB > OP_Ask);           
                        if CR_Target < 1
                            DEV_LB = CR_Target * sum(OP_LB(Index_Ask_W_H) - OP_Ask(Index_Ask_W_H)) + ...
                                     (1 - CR_Target) * sum(OP_Ask(Index_Ask_W_L) - OP_LB(Index_Ask_W_L));
                        else
                            DEV_LB = sum(OP_Ask(Index_Ask_W_L) - OP_LB(Index_Ask_W_L));                        
                        end
                        CR_LB = sum(Index_Ask_W_L) / length(K);
                        clear Index_Ask_W_L Index_Ask_W_H                       
                    
                        Record_CR_Real(c) = CR_LB;
                        Record_DEV_Real(c) = DEV_LB;
                        clear OP_LB CR_LB DEV_LB                       
                    end
                    
                    [AllDEV, ~, Index_AllDEV] = unique(Record_DEV_Real);    
                    if sum(Index_AllDEV==1) > 1
                        break
                    else
                        [~, Index_Min] = min(Record_DEV_Real);
                        Index_Real = Index_Real(Index_Min);                           
                    end
                    clear AllDEV Index_AllDEV

                    OPT_TAB(r, :, d) = [CR_Target Record(Index_Real, 1:2) CR_Real Record_CR_Real(Index_Min) Record_DEV_Real(Index_Min)];
                    clear Index_Min Record_CR_Real Record_DEV_Real
                    
                    OP_LB = Record_OP_LB(Index_Real, :)';   
                    OP_LB = max(OP_LB, 0);                                 
                    OPT_OP_LB(r, d) = cellstr(char(num2str(OP_LB')));
                    OP_UB = Record_OP_UB(Index_Real, :)';   
                    OP_UB = max(OP_UB, 0);                                 
                    OPT_OP_UB(r, d) = cellstr(char(num2str(OP_UB')));
                    clear OP_LB IV_LB OP_UB IV_UB                      
                end
                clear Index_EU
            else
            end
            clear Index_Real
        else
        end

        % Based on Upper Bound
        if strcmp(BV_Target, 'UB')==1
            Index_Real = find(Record(:, end - 1)==CR_Real);                
           if length(Index_Real) > 0
                Index_EU = find((Record(Index_Real, 1)==1) & ...
                                (Record(Index_Real, 2)==1));     
                if length(Index_EU) > 0
                    Index_Real = Index_Real(Index_EU);                           
                    OPT_TAB(r, :, d) = [CR_Target Record(Index_Real, 1:2) CR_Real nan nan];
                    OPT_OP_LB(r, d) = EU_OP_LB(1, d);
                    OPT_OP_UB(r, d) =  EU_OP_UB(1, d);
                else
                    Record_CR_Real = nan(size(Index_Real));                
                    Record_DEV_Real = nan(size(Index_Real));                   
                    for c = 1:length(Index_Real)            
                        OP_UB = Record_OP_UB(Index_Real(c), :)';   
                        OP_UB = max(OP_UB, 0);                             

                        Index_Bid_W_L = logical(OP_UB >= OP_Bid); 
                        Index_Bid_W_H = logical(OP_UB < OP_Bid);           
                        if CR_Target < 1
                            DEV_UB = CR_Target * sum(OP_Bid(Index_Bid_W_H) - OP_UB(Index_Bid_W_H)) + ...
                                     (1 - CR_Target) * sum(OP_UB(Index_Bid_W_L) - OP_Bid(Index_Bid_W_L));
                        else
                            DEV_UB = sum(OP_UB(Index_Bid_W_L) - OP_Bid(Index_Bid_W_L));                        
                        end
                        CR_UB = sum(Index_Bid_W_L) / length(K);
                        clear Index_Bid_W_L Index_Bid_W_H                      

                        Record_CR_Real(c) = CR_UB;
                        Record_DEV_Real(c) = DEV_UB;
                        clear OP_UB CR_UB DEV_UB                       
                    end
                    
                    [AllDEV, ~, Index_AllDEV] = unique(Record_DEV_Real);    
                    if sum(Index_AllDEV==1) > 1
                        break
                    else
                        [~, Index_Min] = min(Record_DEV_Real);
                        Index_Real = Index_Real(Index_Min);                           
                    end
                    clear AllDEV Index_AllDEV

                    OPT_TAB(r, :, d) = [CR_Target Record(Index_Real, 1:2) CR_Real Record_CR_Real(Index_Min) Record_DEV_Real(Index_Min)];
                    clear Index_Min Record_CR_Real Record_DEV_Real
                    
                    OP_LB = Record_OP_LB(Index_Real, :)';   
                    OP_LB = max(OP_LB, 0);                                  
                    OPT_OP_LB(r, d) = cellstr(char(num2str(OP_LB')));
                    OP_UB = Record_OP_UB(Index_Real, :)';   
                    OP_UB = max(OP_UB, 0);                                 
                    OPT_OP_UB(r, d) = cellstr(char(num2str(OP_UB')));
                    clear OP_LB IV_LB OP_UB IV_UB                      
                end
                clear Index_EU
            else
            end
            clear Index_Real
        else
        end
        
        % Based on Both Bounds
        if strcmp(BV_Target, 'LBUB')==1
            Index_Real = find(Record(:, end)==CR_Real);                    
            if length(Index_Real) > 0
                Index_EU = find((Record(Index_Real, 1)==1) & ...
                                (Record(Index_Real, 2)==1));     
                if length(Index_EU) > 0
                    Index_Real = Index_Real(Index_EU);                          
                    OPT_TAB(r, :, d) = [CR_Target Record(Index_Real, 1:2) CR_Real nan nan];
                    OPT_OP_LB(r, d) = EU_OP_LB(1, d);
                    OPT_OP_UB(r, d) =  EU_OP_UB(1, d);
                else
                    Record_CR_Real = nan(size(Index_Real));                
                    Record_DEV_Real = nan(size(Index_Real));               
                    for c = 1:length(Index_Real)   
                        OP_LB = Record_OP_LB(Index_Real(c), :)';   
                        OP_LB = max(OP_LB, 0);                             

                        OP_UB = Record_OP_UB(Index_Real(c), :)';   
                        OP_UB = max(OP_UB, 0);                             
                        
                        Index_Ask_W_L = logical(OP_LB <= OP_Ask);
                        Index_Ask_W_H = logical(OP_LB > OP_Ask);           
                        Index_Bid_W_L = logical(OP_UB >= OP_Bid); 
                        Index_Bid_W_H = logical(OP_UB < OP_Bid);           
                        if CR_Target < 1                                                     
                            DEV_LBUB = CR_Target * (sum(OP_LB(Index_Ask_W_H) - OP_Ask(Index_Ask_W_H)) + ...
                                                    sum(OP_Bid(Index_Bid_W_H) - OP_UB(Index_Bid_W_H))) + ...
                                       (1 - CR_Target) * (sum(OP_Ask(Index_Ask_W_L) - OP_LB(Index_Ask_W_L)) + ...
                                                          sum(OP_UB(Index_Bid_W_L) - OP_Bid(Index_Bid_W_L)));                            
                        else                       
                            DEV_LBUB = sum(OP_Ask(Index_Ask_W_L) - OP_LB(Index_Ask_W_L)) + ...
                                       sum(OP_UB(Index_Bid_W_L) - OP_Bid(Index_Bid_W_L));
                        end
                        CR_LBUB = (sum(Index_Ask_W_L) + sum(Index_Bid_W_L)) / (2 * length(K));
                        clear Index_Ask_W_L Index_Ask_W_H Index_Bid_W_L Index_Bid_W_H                    

                        Record_CR_Real(c) = CR_LBUB;
                        Record_DEV_Real(c) = DEV_LBUB;
                        clear OP_LB OP_UB CR_LBUB DEV_LBUB                       
                    end
                    
                    [AllDEV, ~, Index_AllDEV] = unique(Record_DEV_Real);    
                    if sum(Index_AllDEV==1) > 1
                        break
                    else
                        [~, Index_Min] = min(Record_DEV_Real);
                        Index_Real = Index_Real(Index_Min);                           
                    end
                    clear AllDEV Index_AllDEV

                    OPT_TAB(r, :, d) = [CR_Target Record(Index_Real, 1:2) CR_Real Record_CR_Real(Index_Min) Record_DEV_Real(Index_Min)];
                    clear Index_Min Record_CR_Real Record_DEV_Real
                    
                    OP_LB = Record_OP_LB(Index_Real, :)';   
                    OP_LB = max(OP_LB, 0);                                 
                    OPT_OP_LB(r, d) = cellstr(char(num2str(OP_LB')));
                    OP_UB = Record_OP_UB(Index_Real, :)';   
                    OP_UB = max(OP_UB, 0);                                 
                    OPT_OP_UB(r, d) = cellstr(char(num2str(OP_UB')));
                    clear OP_LB IV_LB OP_UB IV_UB                      
                end
                clear Index_EU
            else
            end
            clear Index_Real
        else
        end
        
        clear CR_Target CR_Real
    end
    clear Record Record_OP_LB Record_OP_UB OP_Bid OP_Ask IV_Bid IV_Ask
    clear K KS CPFlag S0 S0_ADJ TTM RF DY Type
    clear Target_Date Target_TTM   
end
Index_OPT = 1 + [1 2];
Index_CR = 1;
clear d r c

%% Output
if strcmp(BV_Target, 'LB')==1
    FileName = ['Summary_OPT_AllCR_AllPeriod_OP_LB'];
else
end

if strcmp(BV_Target, 'UB')==1
    FileName = ['Summary_OPT_AllCR_AllPeriod_OP_UB'];
else
end

if strcmp(BV_Target, 'LBUB')==1
    FileName = ['Summary_OPT_AllCR_AllPeriod_OP_LBUB'];
else
end

save(FileName)
clear FileName 
